La estructura genética espacial de una metapoblación se caracteriza por el número de subpoblaciones en su interior, las diferencias en frecuencias alélicas entre subpoblaciones y el grado de aislamiento genético de las subpoblaciones. Se dice que hay estructura en las poblaciones, cuando hay diferencias en frecuencias alélicas entre las sub-poblaciones de la meta-población.
Para determinar la estructura genética, se utilizan diferentes estimativas. Entre estas, el índice de fijación FST y GST que miden la diferenciación entre poblaciones, la D de Jost una medida de distancia genética y el Análisis Molecular de Varianza que divide la varianza molecular en sus componentes entre grupos de poblaciones, entre las poblaciones dentro de los grupos y dentro de las poblaciones. Además, existen métodos para visualizar gráficamente y estudiar la relación entre las poblaciones en estudio.
Objetivo general
El objetivo de este laboratorio es conocer las funciones de las librerías en R que permiten el cálculo de las estimativas de estructura genética.
Metodología
Para este laboratorio se utilizará la base de datos magen.csv que consiste en el análisis de 224 muestras de maíz criollo proveniente de las regiones Brunca y Chorotega de Costa Rica, analizadas con 20 microsatélites (SSR). La base de datos está estructurada de manera que se tienen sitios de muestreo por cada región y diferentes colores de grano por sitio de muestreo.
Fig 1: Mapa de Costa Rica con las regiones Chorotega (al norte) y Brunca (al sur) resaltadas en gris. En las regiones se muestran de color los sitios de colecta. Dentro de la Región Chorotega cada círculo corresponde a una accesión en cada sitio de colecta; en la Región Brunca, las accesiones se representan con cuadrados
Para iniciar, deben establecer el directorio de trabajo y cargar las librerías necesarias para hacer la práctica. Se importan los datos a R y luego se generan los estratos, similar a la práctica pasada.
#Recuerden establecer el directorio de trabajo donde están sus datossetwd("/Users/asernas/Desktop/2024/ASISTENCIAS/Genética\ de\ poblaciones\ 2024-1/LAB2-Estructura_F-IBE-IBD/")library(poppr) library(hierfstat) library(pegas) library(tidyr) library(mmod) library(data.table)library(MASS)library(ggplot2)data.maiz <-read.csv("maiz.csv", header =TRUE, sep =";")alelos.maiz <-alleles2loci(data.maiz, ploidy =2, rownames =1, population =2)magen <-loci2genind(alelos.maiz, ploidy =2, na.alleles ="NA")estratos<-read.csv("MaGene_Strata.csv",sep=";",header =TRUE)otros <-list(estratos$REGION, estratos$SITIO, estratos$COLOR)names(otros)<-c("region", "sitio", "color") #cambiar nombresother(magen)<- otros #asignar el objeto al objeto genind magen strata(magen)<-data.frame(magen@other[1:3]) #las columnas 1 a la 3 contienen los estratos de interés debe ser un data.framenameStrata(magen) <-~region/sitio/color magen
/// GENIND OBJECT /////////
// 224 individuals; 20 loci; 351 alleles; size: 396.3 Kb
// Basic content
@tab: 224 x 351 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 4-37)
@loc.fac: locus factor for the 351 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: adegenet::df2genind(X = as.matrix(x[, attr(x, "locicol"), drop = FALSE]),
sep = "/", pop = pop, ploidy = ploidy)
// Optional content
@pop: population of each individual (group size range: 81-143)
@strata: a data frame with 3 columns ( region, sitio, color )
@other: a list containing: region sitio color
1. Estadísticas F
Este grupo de estimativas, surgen a partir de un análisis de diversidad y diferenciación genética de manera jerárquica, donde los individuos (I) se encuentran dentro de sub-poblaciones (S) y las sub-poblaciones forman la población total (T). Las estadísticas F hacen referencia al cálculo de coeficientes de endogamia en estos niveles y entre los mismos. En el cuadro 1 se muestran las tres estadísticas F definidas por Sewall Wright.
Cuadro 1. Descripción y fórmula de las estadísticas F.
Coeficiente
Cálculo
Descripción
FIS
F_{IS} = \frac{H_S - H_I}{H_S}
Coeficiente de endogamia. Mide la reducción en la heterocigosidad de los individuos dentro de las sub-poblaciones que se debe a apareamiento no aleatorio. Va de -1 a 1, donde FIS = -1 exceso de heterocigotos en las sub-poblaciones, FIS = 1, población endógama compuestas por homocigotos.
FST
F_{ST} = \frac{H_T - H_S}{H_T}
Índice de fijación. Mide el déficit de heterocigotas en la metapoblación debido a la pérdida de diversidad genética por deriva, en las sub-poblaciones por el efecto Wahlund. Va de 0 a 1. FST = 0, no hay diferencia en frecuencias alélicas entre sub-poblaciones, no hay efecto Wahlund, no hay estructura. FST = 1, maxima estructura, todas las subpoblaciones están fijadas con alelos distintos.
FIT
F_{IT} = \frac{H_T - H_I}{H_T}
Determina la reducción en heterocigosidad de un individuo relativo a la metapoblación. No es común analizarlo.
Estas estadísticas F propuestas por Sewall Wright fueron originalmente propuestas para un sistema de un único locus con dos alelos. Aunque se ha extendido para ser usado en loci con más de un alelo, no es un estadístico multialélico. Por otro lado, el valor que se reporta de FST cuando se trabaja con más de un locus es el promedio de los distintos loci. Este acercamiento no es del todo correcto, por lo que se han propuesto otros estadísticos multilocus y multialélicos, de los cuales se hablará más adelante.
La librería pegas, permite el cálculo de las tres estimativas por locus. Se debe utilizar el comando Fst(x), donde “x” es un objeto de clase loci. Para esto, es necesario convertir el objeto genind a uno clase loci con la función as.loci(x). Se puede hacer dentro de la misma función Fst(x). Es importante recalcar, que la población o el estrato en que se desea trabajar debe ser definido antes de realizar cualquier cálculo. En el siguiente código, analizaremos la estructura entre sitios dentro de regiones. Es decir, se calculará el valor de FST entre diferentes sitios.
#con pegas por loci. #Medir estructura entre sitios setPop(magen) <-~region/sitio fma<-Fst(as.loci(magen)) fma
Se recomienda asignarlo a un objeto para luego promediar los loci para cada estadístico F. Debido a que el resultado es una matriz con cada loci en filas y las estadísticas F en la columna, se puede utilizar el comando colMeans(x) para obtener los promedios por columna. En este caso al poner que lo haga por sitio según región, el Fst es comparación entre sitios.
colMeans(fma) #promedios por columna
Fit Fst Fis
0.3605690 0.0597368 0.3206533
Este resultado de FST nos indica que el déficit de heterocigotos en la metapoblación causado por el efecto Wahlund es de un 5.97%. Otra forma de interpretarlo es que la variación en frecuencias alélicas entre sitios es de aproximadamente el 5.97%.
El FST se puede obtener de manera pareada (en inglés pairwise), lo que permite estimar la diferencia en frecuencias alélicas entre todos los pares de poblaciones definidas. En este caso, como lo definimos anteriormente vamos a comparar todos los pares de sitios. Se utiliza el comando pairwise.fst.dosage(x, pop=). Se le debe indicar el objeto genind, la población (estrato) sobre la cual se va a trabajar y el tipo de resultado, que se recomienda que sea una matriz de distancias (dist). Aparte del comando pairwise.fst(x), se puede usar el comando pairwise_Gst_Nei(x), que funciona de manera similar al previamente calculado pero utiliza el Gst de Nei, el cual se explicará más adelante.
pf <-pairwise.fst.dosage(magen, pop = magen$strata$sitio) #lo hace por sitiopf
La matriz que se produce es la comparación entre todos los pares de sitios. Hay 10 sitios: “BORUCA”, “LAUREL”, “CONCEPCION”, “VALENCIA”, “STA_CRUZ”, “LIBERIA”, “CANAS”, “CRUZ_LB”, “CRUZ_SC” y “HOJANCHA” . La matriz tiene 10 filas y 10 columnas, sin embargo sólo se muestra la mitad ya que es simétrica. El primer valor de la matriz FST = 0.09316064, es el valor de Fst que compara “BORUCA” con “LAUREL”. Es decir, entre estos dos “sitios” hay un 9.31% de déficit de heterocigotos por efecto Wahlund. Se puede ver que el valor más alto es entre “LAUREL” y “CAÑAS” con un FST = 0.21053541. Estas matrices de comparaciones pareadas, se pueden utilizar para crear árboles o dendrogramas de similitud.
2. GST
El GST es un estadístico propuesto por Nei, que a diferencia del FST, fue desarrollado para marcadores multialélicos. Este es equivalente al FST y se utiliza para medir la diferenciación genética entre poblaciones en relación con la diversidad genética total (al igual que un FST). El comando Gst_Nei(x) de la librería mmod, calcula la estimativa, para cada loci y de manera global, a partir de un objeto genind. GST = 0 indica que no hay diferencia en frecuencias alélicas entre las poblaciones, mientras que GST = 1 indica que las poblaciones están fijadas para diferentes alelos. En muchos artículos, cuando se reporta un FST multilocus, usualmente se calculó con un GST.
#El valor varía según la estratificación indicada con setPop, en este caso está por ~region/sitioGst_Nei(magen)
El comando produce un GST global de GST = 0.076. Este valor indica que el efecto Wahlund es del 7.61%. Dicho de otra forma, que la varianza en frecuencias alélicas entre sitios es de 7.61% de la variación total. Se nota que este estimativo es un poco mayor que el FST calculado anteriormente.
También se puede calcular el G’ST o el GST de Hedrick, utilizando el comando GstHedrick(x). A diferencia del GST de Nei, el GST de Hedrick corrige el estimativo de Nei dividiendolo por el valor del GST máximo teórico basado en la heterocigosidad en ese marcador. De esta manera, el G’ST de Hedrick nos indica la estructura genética como un porcentaje de la estructura máxima posible. Esto permite que se compare mejor entre estudios.
Es importante identificar si los valores de GST son significativos, es decir, probar si el valor de GST es estadísticamente diferente de cero. La función chao_bootstrap(x) implementa un muestreo tipo bootstrap al archivo genind, donde se seleccionan muestras al azar (con remplazo) dentro de cada subpoblación. Este muestreo se realiza acorde al tamaño de la subpoblación. La función summarise_bootstrap(x) aplica las estadísticas de diferenciación a la lista de archivos genind generadas con la función anterior. Es importante reportar, ya sea en materiales y métodos o junto con el valor de probabilidad obtenido, cuántas repeticiones de la simulación se realizaron; la función summarise_bootstrap(x) utiliza por defecto 1000 repeticiones.
Primero, se debe establecer el estrato en el que se quiere realizar la prueba. Luego, se usa la función chao_bootstrap(x). Por default, la cantidad de repeticiones es 1000 y puede ser modificado con el argumento nrep=. Para este laboratorio, por razones de tiempo realizaremos 100 replicas, pero lo comun es 1000 o más. Con la función summarise_bootstrap(x), se debe indicar la lista de archivos genind y el estadístico que se quiere realizar, los cuales pueden ser: D_jost, GST_Hedrick y GST_Nei. Si el intervalo de confianza no incluye el 0, se puede afirmar que se cuenta con estructura significativa.
Estimates for each locus
Locus Mean 95% CI
bnlg1046 0.0703 (0.051-0.090)
bnlg1191 0.0821 (0.064-0.100)
bnlg1194 0.0580 (0.045-0.071)
bnlg1265 0.0535 (0.040-0.067)
bnlg1449 0.1796 (0.151-0.208)
bnlg1523 0.0722 (0.057-0.087)
bnlg1740 0.0598 (0.040-0.080)
bnlg1890 0.0638 (0.047-0.081)
bnlg1917 0.0823 (0.062-0.103)
bnlg2305 0.0711 (0.051-0.091)
phi015 0.0447 (0.010-0.079)
phi031 0.1109 (0.074-0.148)
phi064 0.0581 (0.037-0.079)
phi072 0.0620 (0.041-0.083)
phi078 0.1155 (0.084-0.147)
phi089 0.1498 (0.096-0.204)
phi093 0.0676 (0.046-0.090)
phi109188 0.0400 (0.018-0.063)
phi116 0.0974 (0.062-0.133)
phi96100 0.0131 (-0.012-0.038)
Global Estimate based on average heterozygosity
0.0761 (0.070-0.082)
Este resultado indica el valor de intervalo de confianza para cada locus y el valor final para para la estadística global. Este resultado indica que existe estructura genética significativa a nivel de sitios, con un valor de GST = 0.0761, y un límite de confianza al 95% de CI: [0.070-0.082] con 100 repeticiones. Como el CI no incluye el cero, se identifica que la estructura genética es significativa con un nivel de confianza del 95%.
4. D Jost
La D de Jost es una medida de distancia genética, que calcula la fracción de variación alélica entre las poblaciones. Se puede estimar con el comando D_Jost(x), donde x es un objeto genind. Como resultado se obtiene una lista con la D para cada loci, así como dos estimaciones globales: $Global.het usa los promedios de Hs y Ht en todos los loci mientras que $global.harm_mean toma la media armónica de todos los loci. Ambos valores se utilizan para el cálculo de la D.
Debido a que el objeto $per.locus de la lista DJ es un valor numérico, se puede obtener el promedio de la D de Jost utilizando la función mean().
mean(DJ$per.locus) #promedio de D Jost
[1] 0.3385399
Con el valor obtenido, se puede interpretar que existe una distancia genética por D de Jost de 33.9% promedio entre los sitios muestreados.
2. IBD, IBE y NJ
Introducción
Cuando una población cuenta con estructura genética espacial, se tiene el interés de conocer cuáles son los factores ambientales, bióticos o geográficos que pueden afectar la estructuración y el aislamiento de las subpoblaciones. Para ello, se pueden hacer correlaciones entre los valores de diferenciación genética entre las poblaciones y diferentes variables, como por ejemplo la distancia entre las poblaciones, las diferencias en elevación, temperatura, precipitación, u otros aspectos ambientales.
Existen métodos para visualizar gráficamente y estudiar la relación entre las poblaciones en estudio y los factores que puedan afectar la estructuración espacial. Uno de los métodos más utilizados es la construcción de dendrogramas a partir de las diferencias genéticas entre las poblaciones.
Objetivo general
Aprender a calcular correlaciones entre las estimativas de diferenciación genética y variables abióticas, así como aprender a visualizar la distancia genética entre poblaciones.
1. Aislamiento por distancia (IBD)
El aislamiento por distancia (IBD) propone que existe una correlación entre las distancias geográficas y las distancias genéticas. De manera que las poblaciones más cercanas entre sí, son más similares genéticamente. Para realizar la prueba de IBD se calculan las diferencias genéticas pareadas entre todas las poblaciones. De igual forma se estima una distancia entre poblaciones. Esta distancia puede ser medida como una línea entre poblaciones, o la distancia euclidiana entre ellas. Las diferencias geográficas y las genéticas se comparan con la prueba de correlación de Mantel que compara dos matrices. La prueba aleatoriza los datos, y compara el coeficiente de correlación calculado, con la distribución de coeficientes de correlación calculados a partir de datos aleatorizados (i.e., donde no debería haber correlación). Si el valor estimado es mayor al 95% de las simulaciones, se asume que hay una correlación significativa entre la distancia geográfica y la genética. Es decir, entre más separadas estén las poblaciones, mayor será la diferencia en frecuencias alélicas.
La matriz de distancia geográfica se obtiene con el comando dist(x) a partir de una tabla de datos que contiene las coordenadas (x, y), (longitud, latitud) de los sitios de interés. Es importante que el orden de los sitios en la tabla de datos con las coordenadas, sea igual que el orden de los sitios de la matriz de datos genéticos. El comando dist(x) produce distancias euclidianas entre las coordenadas. Estas distancias no toman en cuenta la topografía o la esfericidad del planeta, por lo que para puntos muy distanciados pueden contener errores. En Costa Rica las distancias son tan cortas, que no representan un error.
Lo primero que vamos a hacer es importar una tabla de datos con los sitios y las coordenadas geográficas, las cuales se denominan x y y respectivamente. La opción row.names = 1 le indica a R que la columna número 1, incluye los nombres de las filas y sep = ","` indica que las columnas están separadas por comas.
dgmaiz <-read.csv("datos_dist.csv", row.names=1, header =TRUE, sep =",") #importa tabla de datos head(dgmaiz) #visualizar los primeros datos
x y
BORUCA -83.34250 9.015347
LAUREL -82.91706 8.460722
CONCEPCION -83.48820 9.097384
VALENCIA -83.72333 9.446394
STA_CRUZ -85.59671 10.219330
LIBERIA -85.45480 10.735898
Note que en la primera columna está la longitud en grados. El negativo indica que se encuentra en longitud oeste. La segunda columna incluye la latitud de las poblaciones.
Dgeo <-dist(dgmaiz) #calcular la matriz de distancias geográficas euclidianas Dgeo
La matriz que produce el comando dist (x) es una matriz triangular, donde se muestra la distancia (en grados) entre cada par de poblaciones. Es decir, entre BORUCA y LAUREL hay 0.6990 grados de distancia. La escala de la distancia no importa mucho, mientras tengamos claro qué se está midiendo. En Costa Rica se pueden usar otras unidades de posición como coordenadas Lambert, o UTC; de esta forma, las distancias entre poblaciones estarían dadas en kilómetros.
Para las distancias genéticas, podemos utilizar el comando pairwise.fst(x) o el comando pairwise_Gst_Nei(x), que funciona de manera similar al previamente calculado pero utiliza el Gst de Nei. En este caso usaremos el Gst de Nei.
setPop(magen) <-~sitio #indicar que la estratificación es por sitios Dgen <-pairwise_Gst_Nei(magen)Dgen
La matriz Dgen contiene los valores de GST pareados entre poblaciones. Por ejemplo, entre BORUCA y LAUREL hay un GST = 0.053. La mayor diferencia en frecuencias alélicas se encuentra entre las poblaciones de VALENCIA y CONCEPCIÓN, dónde el GST = 0.0967. Lo que interesa ahora es determinar si la matriz de distancias euclidianas está correlacionada con la matriz de distancias genéticas, esto mediante una correlación de Mantel.
Con el siguiente comando, se realiza la correlación de Mantel. Se puede modificar el número de repeticiones con el argumento nrepet=.
#Realizar la prueba de Mantelibd <-mantel.randtest(Dgeo, Dgen, nrepet=999) ibd
Monte-Carlo test
Call: mantel.randtest(m1 = Dgeo, m2 = Dgen, nrepet = 999)
Observation: 0.1852853
Based on 999 replicates
Simulated p-value: 0.126
Alternative hypothesis: greater
Std.Obs Expectation Variance
1.25344280 0.00800654 0.02000343
El valor de la correlación sería 𝑟𝑀 = 0. 1852 y esta correlación no difiere significativamente de cero: p > 0.05 (El valor de p, cambia cada vez que se corre). Por lo que parece que no existe aislamiento por distancia. De manera gráfica, se puede observar la correlación entre las matrices de distancia genética (y) y geográfica (x).
df<-data.frame(dist_geo =as.vector(Dgeo),dist_gen =as.vector(Dgen))ggplot(df, aes(dist_geo, dist_gen))+geom_point()+geom_smooth(method ="lm", se = F, color ="red", linetype ="dashed")+theme_classic() +theme(axis.text.y =element_text(angle =90, vjust =0.5, hjust =1))
`geom_smooth()` using formula = 'y ~ x'
2. Aislamiento por ambiente (IBE)
El aislamiento por ambiente (IBE por sus siglas en inglés: Isolation by environment) propone que existe un incremento en las diferencias genéticas entre las poblaciones a medida que aumenta la diferencia ambiental entre ambas poblaciones independiente de su distancia geográfica. Estos pueden ser factores abióticos, como elementos climáticos, elevación, temperatura, o bióticos como densidad vegetal, tipo de ecosistema, entre otros. En su mayoría, estos factores se miden de forma continua, lo cual permite una relación más simple entre la diferenciación genética y la ambiental, aunque también pueden usarse factores discretos.
Al igual que con IBD, para estimar IBE se puede realizar una prueba mantel cuando los datos ambientales son continuos. En esta práctica, se utilizará como variable ambiental la elevación de la ubicación de cada población de la base de datos de maíz silvestre. Al igual que con IBD, se debe tener una base de datos (archivo) con la elevación de cada población y estimar la diferencia entre ellas, para luego realizar la prueba mantel al correlacionarse con las diferencias genéticas. Es importante mencionar que la elevación se puede incluir como una variable además de la latitud y la longitud. Es decir, uno puede calcular la distancia euclidiana entre dos puntos utilizando múltiples dimensiones. Algunos autores sugieren utilizar múltiples variables ambientales para cada sitio y mediante un análisis de componentes principales, estimar distancias entre sitios. En el ejemplo de este laboratorio, haremos una comparación simple, dónde se calculará la distancia entre sitios sólo usando la elevación.
setPop(magen) <-~sitioaltitud=read.csv2(file ="datos_altitud.csv") #cargamos la base de datosrownames(altitud)=altitud$sitio #asignamos el nombre de los sitiosDalt=dist(altitud[2]) #calculamos la distancia euclidiana entre las altitudesDgen=pairwise_Gst_Nei(magen) mantel.randtest(Dgen,Dalt)
Monte-Carlo test
Call: mantel.randtest(m1 = Dgen, m2 = Dalt)
Observation: 0.6646986
Based on 999 replicates
Simulated p-value: 0.01
Alternative hypothesis: greater
Std.Obs Expectation Variance
2.28426478 0.01894211 0.07991801
Acá obtenemos un valor de correlación alto, de 𝑟𝑀 = 0. 664, y esta correlación difiere significativamente de cero con 999 repeticiones ( p < 0.05 ). Esto implica que existe aislamiento causado por la elevación, en donde las poblaciones encontradas a mayor elevación serán distintas genéticamente a aquellas a más baja elevación. Podemos visualizar esta relación de la misma forma que IBD.
df<-data.frame(dist_alt =as.vector(Dalt),dist_gen =as.vector(Dgen))ggplot(df, aes(dist_alt, dist_gen))+geom_point()+geom_smooth(method ="lm", se = F, color ="red", linetype ="dashed")+theme_classic() +theme(axis.text.y =element_text(angle =90, vjust =0.5, hjust =1))
`geom_smooth()` using formula = 'y ~ x'
Aunque tenemos una correlación significativa, nuestro muestreo presenta un sesgo de muestreo por poblaciones a baja altitud. Para poder concluir con mayor certeza que existe un efecto de la elevación sobre la diferenciación genética, es recomendable aumentar el muestreo de sitios a elevación intermedia y alta.
Si comparamos los valores de correlación, podemos afirmar que existe un mayor efecto de IBE que de IBD sobre las poblaciones de maíz silvestre muestreadas.
3. Dendrograma
Con las distancias se puede generar dendrogramas por el algoritmo de neighbour-joining (NJ: https://en.wikipedia.org/wiki/Neighbor_joining), que va generando el árbol permutando comparaciones y uniendo las poblaciones que son más similares entre sí. Para esto se utiliza una serie de comandos, entre estos nj(x) que opera sobre una matriz de distancias (x) y realiza un árbol; ladderize() que reorganiza la estructura interna del árbol para obtener un efecto escalonado. Finalmente, el árbol se grafica con el comando plot().
setPop(magen) <-~region/sitio #se establece el estratoarbol.sitio <-pairwise_Gst_Nei(magen)%>%nj() %>%# calcula el árbol njladderize() # organiza las ramas por poblaciónplot(arbol.sitio) #realiza el gráfico
Se recomienda otras formas de visualizar el mismo árbol, debido a que no se están tratando relaciones filogenéticas.
#tipos de árbolplot(arbol.sitio, type="unrooted", font=1, cex=0.5)
En estos dendrogramas, las poblaciones que se encuentran unidas o están más cerca en el árbol, son aquellas que tienen menor diferencias en frecuencias alélicas. Es decir, menos estructura entre ellas.
Se puede realizar el árbol para otros niveles de estratificación, solo se debe indicar en setPop() cuál estrato se está trabajando y asegurarse de que las distancias genéticas sean calculadas según el nuevo estrato.
